Main
Data processing
Select conditions with minimum replicates
filtered.samples.retained <- list()
filtered.samples.retained$conditions <- filtered.samples$conditions %>%
dplyr::filter(n>=params$min.replicates)
filtered.samples.retained$metadata <- filtered.samples$metadata %>%
dplyr::filter(grepl(paste(filtered.samples.retained$conditions$regex,collapse="|"),.$SampleName))
if ("Sex" %in% colnames(filtered.samples.retained$conditions)) {
conditions <-
paste(filtered.samples.retained$conditions$Tissue,
filtered.samples.retained$conditions$Sex, sep="_")
} else {
conditions <- filtered.samples.retained$conditions$Tissue
}
filtered.samples.retained$log2.tmm.cpm.gf <- vector(mode="list", length=length(conditions))
names(filtered.samples.retained$log2.tmm.cpm.gf) <- conditions
jack.condition.retained <- vector(mode="list", length=length(conditions))
names(jack.condition.retained) <- conditions
for (c in 1:length(conditions)) {
if ("Sex" %in% colnames(filtered.samples.retained$conditions)) {
t <- filtered.samples.retained$conditions$Tissue[[c]]
s <- filtered.samples.retained$conditions$Sex[[c]]
ts <- paste(t,s,sep="_")
filtered.samples.retained$log2.tmm.cpm.gf[[ts]] <- filtered.samples$log2.tmm.cpm.gf[[t]][[s]]
jack.condition.retained[[ts]] <- jack.bind %>%
dplyr::filter(Biotype=="protein_coding") %>%
dplyr::filter(paste(.$Tissue, .$Sex, sep="_")==ts)
} else {
ts <- filtered.samples.retained$conditions$Tissue[[c]]
filtered.samples.retained$log2.tmm.cpm.gf[[ts]] <- filtered.samples$log2.tmm.cpm.gf[[ts]]
jack.condition.retained[[ts]] <- jack.bind %>%
dplyr::filter(Biotype=="protein_coding") %>%
dplyr::filter(Tissue==ts)
}
# Filter gene list
filtered.samples.retained$log2.tmm.cpm.gf[[ts]] <- filtered.samples.retained$log2.tmm.cpm.gf[[ts]] %>%
.[rownames(.) %in% jack.condition.retained[[ts]]$GeneID,]
}
Compute z-score per condition
zscore.matrix <- lapply(setNames(conditions, conditions), function(ts) {
zscore(filtered.samples.retained$log2.tmm.cpm.gf[[ts]])
})
Compute bimodality test per condition
bimodality.test <- lapply(setNames(conditions, conditions), function(ts) {
print(ts)
bimodalIndex(zscore.matrix[[ts]]) %>%
tibble::rownames_to_column(var="GeneID")
})
## [1] "brain_F"
## 1 ..........
## 2 ..........
## 3 ..........
## 4 ..........
## 5 ..........
## 6 ..........
## 7 ..........
## 8 ..........
## 9 ..........
## 10 ..........
## 11 ..........
## 12 ..........
## 13 ..........
## 14 .........
## [1] "brain_M"
## 1 ..........
## 2 ..........
## 3 ..........
## 4 ..........
## 5 ..........
## 6 ..........
## 7 ..........
## 8 ..........
## 9 ..........
## 10 ..........
## 11 ..........
## 12 ..........
## 13 ..........
## 14 ......
## [1] "eye_F"
## 1 ..........
## 2 ..........
## 3 ..........
## 4 ..........
## 5 ..........
## 6 ..........
## 7 ..........
## 8 ..........
## 9 ..........
## 10 ..........
## 11 ..........
## 12 ..........
## 13 ..........
## 14 .....
## [1] "eye_M"
## 1 ..........
## 2 ..........
## 3 ..........
## 4 ..........
## 5 ..........
## 6 ..........
## 7 ..........
## 8 ..........
## 9 ..........
## 10 ..........
## 11 ..........
## 12 ..........
## 13 ..........
## 14 .......
## [1] "gonads_F"
## 1 ..........
## 2 ..........
## 3 ..........
## 4 ..........
## 5 ..........
## 6 ..........
## 7 ..........
## 8 ..........
## 9 ..........
## 10 ..........
## 11 ..........
## 12 ..........
## 13 ....
## [1] "intestine_F"
## 1 ..........
## 2 ..........
## 3 ..........
## 4 ..........
## 5 ..........
## 6 ..........
## 7 ..........
## 8 ..........
## 9 ..........
## 10 ..........
## 11 ..........
## 12 ..........
## 13 ....
## [1] "intestine_M"
## 1 ..........
## 2 ..........
## 3 ..........
## 4 ..........
## 5 ..........
## 6 ..........
## 7 ..........
## 8 ..........
## 9 ..........
## 10 ..........
## 11 ..........
## 12 ..........
## 13
## [1] "liver_F"
## 1 ..........
## 2 ..........
## 3 ..........
## 4 ..........
## 5 ..........
## 6 ..........
## 7 ..........
## 8 ..........
## 9 ....
## [1] "pectoral_fin_F"
## 1 ..........
## 2 ..........
## 3 ..........
## 4 ..........
## 5 ..........
## 6 ..........
## 7 ..........
## 8 ..........
## 9 ..........
## 10 ..........
## 11 ..........
## 12 ..........
## 13 .......
## [1] "pectoral_fin_M"
## 1 ..........
## 2 ..........
## 3 ..........
## 4 ..........
## 5 ..........
## 6 ..........
## 7 ..........
## 8 ..........
## 9 ..........
## 10 ..........
## 11 ..........
## 12 ..........
## 13 ........
## [1] "skin_F"
## 1 ..........
## 2 ..........
## 3 ..........
## 4 ..........
## 5 ..........
## 6 ..........
## 7 ..........
## 8 ..........
## 9 ..........
## 10 ..........
## 11 ..........
## 12 ..........
## 13 ..........
## 14
## [1] "skin_M"
## 1 ..........
## 2 ..........
## 3 ..........
## 4 ..........
## 5 ..........
## 6 ..........
## 7 ..........
## 8 ..........
## 9 ..........
## 10 ..........
## 11 ..........
## 12 ..........
## 13 .......
Bin variability ranks
for (ts in conditions) {
bimodality.test[[ts]]$VarRank <-
vapply(bimodality.test[[ts]]$GeneID, FUN.VALUE=numeric(1), FUN=function(g) {
jack.condition.retained[[ts]]$Mean_Local_Rank_Log2CV[jack.condition.retained[[ts]]$GeneID==g]
})
# Group variability ranks into bins of size 0.1
bimodality.test[[ts]]$VarRankBin <- round(bimodality.test[[ts]]$VarRank, 1)
}
Plots
Boxplots - Bimodality index (BI) by variability rank
# Bind data frames for plotting
bimodality.test.bind <- dplyr::bind_rows(bimodality.test, .id="Condition")
Landscape
boxplot_bimodality_by_variability(bimodality.test.bind, title=params$species.name, facet=TRUE)
Portrait
boxplot_bimodality_by_variability(bimodality.test.bind, title=params$species.name, facet=TRUE)
Expression level z-scores by variability rank bin
Question: At which variability rank does the z-score distribution visually appear bimodal?
All
# Will only print last plot to HTML, but works on RStudio
for (ts in conditions) {
plot_density_zscore(zscore.matrix[[ts]], bimodality.test[[ts]], title=ts)
}
Example: Brain
if ("brain" %in% filtered.samples.retained$conditions$Tissue) {
if ("Sex" %in% colnames(filtered.samples.retained$conditions)) {
ts <- "brain_F"
} else { ts <- "brain" }
plot_density_zscore(zscore.matrix[[ts]], bimodality.test[[ts]], title=ts)
}
Example: Gonads
if ("gonads" %in% filtered.samples.retained$condtions$Tissue) {
if ("Sex" %in% colnames(filtered.samples.retained$conditions)) {
ts <- "gonads_F"
} else { ts <- "gonads" }
plot_density_zscore(zscore.matrix[[ts]], bimodality.test[[ts]], title=ts)
}